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ABSTRACT 

The nonlinear force-free field (NLFFF) model is often used to describe the solar coronal magnetic field, however 
a series of earlier studies revealed difficulties in the numerical solution of the model in application to photospheric 
boundary data. We investigate the sensitivity of the modeling to the spatial resolution of the boundary data, by applying 
multiple codes that numerically solve the NLFFF model to a sequence of vector magnetogram data at different resolu¬ 
tions, prepared from a single Hinode/SOT-SP scan of NOAA Active Region 10978 on 2007 December 13. We analyze 
the resulting energies and relative magnetic helicities, employ a Flelmholtz decomposition to characterize divergence 
errors, and quantify changes made by the codes to the vector magnetogram boundary data in order to be compatible with 
the force-free model. This study shows that NLFFF modeling results depend quantitatively on the spatial resolution of 
the input boundary data, and that using more highly resolved boundary data yields more self-consistent results. The free 
energies of the resulting solutions generally trend higher with increasing resolution, while relative magnetic helicity 
values vary significantly between resolutions for all methods. All methods require changing the horizontal components, 
and for some methods also the vertical components, of the vector magnetogram boundary field in excess of nominal 
uncertainties in the data. The solutions produced by the various methods are significantly different at each resolution 
level. We continue to recommend verifying agreement between the modeled field lines and corresponding coronal loop 
images before any NLFFF model is used in a scientific setting. 


Subject headings: Sun: corona — Sun: magnetic fields 

1. Introduction 

The solar coronal magnetic field produces solar activity, 
including extremely energetic solar flares and coronal mass 
ejections. There is considerable interest in accurate modeling 
of magnetic fields in and around active regions on the Sun, the 
locations of the most intense coronal fields and the drivers of 
flares and many mass ejections, with the aim of better under- 
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Standing the physics underlying magnetic energy release, and 
improving the ability to predict space weather storms caused 
by large events. 


A popular model for the coronal magnetic held B is the 
nonlinear force-free held (NLFFF) model (see the reviews by 
Wiegelmann & Sakurai|[2nT2 and Regnier|[20l3| l, which as¬ 


sumes a static configuration with a zero Lorentz force, 


JxB = 0, (I) 

where J = V X B is the electric current density, together 
with the solenoidal condition. 


V • B = 0. (2) 

The model current density is everywhere parallel to the mag¬ 
netic held, and Equation ([T]) is often written 

V X B = aB, (3) 

where a is the force-free parameter. The proper boundary 
conditions on the model are the specification of the normal 
component of the held Bn over the bounding surfaces of the 
solution domain, together with the normal component of the 
electric current density J„ (or, alternatively, a) over one po¬ 
larity of the held in the boundary (i.e., either over the region 
where Bn > 0 or over the region where Bn < 0) (|Grad &| 
Rubin||T958] l. The current density is only required over one 
polarity of the held because current density streamlines fol¬ 
low magnetic held lines, so values of J„ at one polarity are 
“mapped” to the other polarity by the geometry of the held 













lines of the solution. Stated another way, taking the diver¬ 
gence of Equation ([^ and applying the solenoidal condition 
yields the general property 


B • Va = 0, 


(4) 


which indicates that a is invariant along magnetic field lines 
(though a may be different from line to line). In the follow¬ 
ing, we consider the problem in a half space (z > 0), with 
the plane z = 0 representing the base of the (assumed force- 
free) corona. This geometry neglects solar curvature, which 
is appropriate for the local modeling done in this study. 

Solar observations provide a set of boundary data for ap¬ 
plication of the model. Spectro-polarimetric measurements of 
magnetically sensitive lines ( |del Toro Iniesta|2003| l are used 
to constmct vector (three-component) magnetogram data that 
are presumed to be located on a planar surface representing 
the top of the photosphere. The vertical current density 
may be estimated from such vector magnetogram data using 


f^oJz\z=0 ~ 


dB„ dB^, 


dx 


dy 


(5) 


z=0 


provided the 180-degree ambiguity in the direction of the held 
transverse to the line of sight is resolved ( Metcalf et al.|2006 
|Leka et al.|2009| l. 

Vector magnetogram data are now routinely produced by 
ground based instruments, and more recently have been pro¬ 
vided by two space-based vector magnetographs: the So¬ 
lar Optical Telescope/Spectro-Polarimeter (SOT-SP) on the 
Hinode satellite ( |Tsuneta et al.|[2008| ), and the Helioseis- 
mic and Magnetic Imager on the Solar Dynamics Observa¬ 
tory (SDO/HMI; |Schou et al.||20I2} . The use of such vec¬ 
tor magnetogram data for NLFFF modeling builds on earlier 
work done with data from the Haleakala Stokes Polarimeter 
(HSP), the Imaging Vector Magnetograph (IVM), the Solar 
Flare Telescope (SFT), and the Advanced Stokes Polarime¬ 
ter (ASP) in the previous decades (e.g., Mikic & McClymont| 
1994 Roumeliotis I996[ Thalmann & Wiegelmann 20081. 


The earliest vector magnetogram data possessed spatial res¬ 
olutions of multiple arc seconds, whereas the Hinode/SOT-SP 
magnetogram data have a resolution as high as 0'.'32. 

In practice, however, additional modeling assumptions are 
needed in order to use vector magnetogram data with NFFFF 
modeling. One problem is that the boundary conditions on 
(or a) are inconsistent with the NFFFF model over the two 
polarities of B^ ( Molodenskii| 1 96^ Aly|1984 19891. Addi¬ 
tionally, boundary data are not available at the top and side 
surfaces of the three-dimensional solution domain. Further¬ 
more, the vector magnetogram measurements contain uncer¬ 
tainties. The polarization measurements are subject to obser¬ 
vational uncertainty, and the process of determining magnetic 
held values by inverting the Stokes polarization spectra in¬ 
volves making various assumptions about the radiative trans¬ 
port of polarized radiation through the magnetized solar at¬ 
mosphere. 

Given the scientihc importance of determining the free en¬ 
ergy in the solar coronal magnetic held, coupled with the 


recent increase in availability of vector data over the past 
decade, a sequence of yearly workshops was organized be¬ 
tween 2004 and 2009 in an effort to characterize and improve 
NFFFF modeling. The workshops demonstrated that NFFFF 
methods work for analytic test cases and for synthetic, solar¬ 
like test data, but encounter specihc problems in application to 
photospheric vector magnetogram data ( SchrijveretalJ20^ 


[Metcalf et al.|200^|Schrij ver et al.|2008| De Rosa et al.|200^ . 

There are a number of different methods of solution of the 
NFFFF model in use, which were found to produce signih- 
cantly different results. Discrepancies include the locations 
and magnitudes of currents within the solution volume, and 
the total magnetic energy in the solution domain. Issues iden- 
tihed during these studies include the inconsistency of photo- 
spheric vector magnetogram data with the NFFFF model, the 
possibility that high spatial-resolution data are needed to ac¬ 
count for small-scale currents, the limited held of view of the 
data, and the lack of account of the substantial uncertainties 
in the boundary held values. 

More recently, additional NFFFF modeling workshop^ 
were held to address these issues. In this article, we follow up 
on one such issue and characterize the inhuence of the spatial 
resolution of vector magnetogram data on the results of the 
modeling. Examining this issue at this time is motivated by 
the increasing availability of vector data from various space- 
and ground-based instrumentation, all of which provide vec¬ 
tor data at different spatial resolutions. For instance, the Na¬ 
tional Solar Observatory’s Synoptic Optical Fong-term Inves¬ 
tigations of the Sun Vector SpectroMagnetograph instrument 
(SOFIS/VSM) provides full disk data with a spatial sampling 
of 1"1 ( jHenney et al.||2006||200^ , while SDO/H MI full-disk 
vector magnetograms have a pixel size of Of.'S (jHoekserna 


et al.|2014| , and “normal-map” and “fast-map” Hinode/SOT- 


SP data are sampled at 0'.'16 and (y/32 respectively (Fites et al. 
|20T3] |. 


To assess the sensitivity of the results of NFFFF model¬ 
ing to variations in spatial resolution, we use data for NOAA 
Active Region (AR) 10978. Vector magnetograms are con- 
stracted from a Hinode/SOT-SP normal-map scan of this re¬ 
gion, with the spatial resolution of the data artihcially de¬ 
graded by a sequence of binning factors. The methodology 
is to perform inversions on rebinned polarization spectra (as 
opposed to simply rebinning the resulting vector data inverted 
from spectra at the native Hinode/SOT-SP resolution) in order 
to approximate observations of AR 10978 by instraments hav¬ 
ing different spatial resolutions. In common with the earlier 
workshop studies, we apply a number of different methods 
of solution of the NFFFF model in order to also gauge the 
dependence of the results on the solution method. 


The effects of spatial resolution on coronal held modeling 
have been discussed in several earlier studies. |Parker| ( ri996| ) 
claimed that the concentration of photospheric magnetic helds 


^ These workshops were hosted and in part supported hy the International 
Space Science Institute (ISSI) in Bern, Switzerland. The first workshop was 
held from 2013 January 29 to February 1 and the second 2014 January 13-16. 
Some online meeting materials are available at http : / /www. issibern. 
ch/teams/solarcorona 
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into unresolved fibrils renders currents inferred from vector 
magnetograms meaningless, but McClymont et al. (1997]) ar¬ 
gued in response that it is only necessary to resolve the large- 
scale twist in the field to correctly infer the current. Semel &| 
Skumanich (1998 i presented a method for inferring | J^j (for 


observations at disk center) that is independent of the ambi¬ 
guity resolution, and argued for the reality of electric currents 
obtained using this method applied to ASP data. More re¬ 
cently, |Leka et ah] ( |2009| l investigated the influence of noise 
and spatial resolution on the ambiguity resolution step in the 
treatment of vector magnetogram data (see also [ GeorgouE^ 
|2012||Leka et al.|2012[ and |Crouch|2013| l. A hare-and-hounds 
exercise was run on an analytic test case, and it was shown 
that failure to resolve observed structures due to low resolu¬ 
tion leads to serious errors in ambiguity resolution. 

|Leka & Barnes | ( |2012| l also looked at the results of arti- 
hcial degradation of Hinode/SOT-SP data on vector magnetic 
field values, including degrading observed Stokes spectra, and 
using full-resolution spectra but degrading the held values de¬ 
rived from the spectra. They examined, in particular, the ef¬ 
fect of these steps on the statistics of the derived helds. They 
found that degraded resolution data can exhibit increased av¬ 
erage Hux densities, lower total hux, and held vectors shifted 
towards the line of sight. The distribution of was found to 


be very sensitive to the spatial resolution. Recently Thalmann 


et al. (2013|l compared NLFFF reconstructions for an active 


region based on Hinode/SOT-SP and SDO/HMl data, includ¬ 
ing calculations for the Hinode data at original resolution and 
rebinned to match the HMI data. They found similar results 
for the different-resolution data from Hinode, but signihcantly 
different results for data between the two instruments, in par¬ 
ticular, e.g., differences in magnetic connectivity. We note 
that Thalmann et al. ( |2013 i rebinned vector magnetogram 
held values, rather than rebinning the original Hinode/SOT- 
SP spectral data as is done in the experiments presented here. 
Rebinning the spectra was shown in Leka & Barnes ( 2012| l to 
produce different vector magnetogram values than rebinning 
the held values. These earlier studies motivate the present in¬ 
vestigation of the inhuence of resolution on NLFFF modeling. 


The structure of the paper is as follows. In Section]^ the 
data used and the different methods are presented, with Sec¬ 
tion | 2 ^ describing the preparation of the different resolution 
vector magnetograms, and Section |2.2| explaining the meth¬ 
ods of solution of the NLFFF model, and their treatment of 
boundary conditions. In Section]^ we present a quantitative 
analysis of the energy and relative magnetic helicity of the 
resulting solution helds (Section [3T]), including assessment 
of non-solenoidal held errors, and an analysis of how the so¬ 
lution methods alter the vector magnetogram boundary data 
(Section jr^. In Section]^ results are discussed and conclu¬ 
sions drawn. 


2. Data and Methods 

2.1. Vector Magnetogram Data for Active Region 10978 

An ideal active region for this study would be one that is 
hux-balanced and is mostly isolated, shows evidence of non¬ 


potentiality, is located near disk center, is small enough for the 
full active region to be observed with high-resolution spec- 
tropolarimetry, and has accompanying high-resolution coro¬ 
nal imagery suitable for comparisons between coronal loops 
and the held lines from the resulting extrapolations. Unfortu¬ 
nately, no such target satisfying all of these criteria was found 
from archive data. As a result, we prioritized the require¬ 
ments for the nonpotentiality, location, isolation, and size for 
the magnetic held, and present in Sectionj^performance met¬ 
rics that do not rely on coronal loop data. Although some 
data from both the X-Ray Telescope (XRT) on Hinode and the 
Transition Region and Coronal Explorer (TRACE) instrument 
were available, there were an insufficient number of distinct 
and discernible loops for useful analysis. 


Based on the above specihcations, NOAA AR 10978 was 
selected from the archive of Hinode/SOT-SP observations. 
The SOT-SP instrument returns full Stokes spectra for the 
Ee I doublet at approximately 630.2 nm with high spectral 
sampling (2.15 pm). Polarization spectra for AR 10978 were 
obtained in normal-map mode on 2007 December 13, 12:18- 
13:41 UTE] 

over a held of view of approximately 164"xl64" (on a 
1024x1024 grid), with high spatial sampling (0'.'16) both 
along the slit and in the direction of the scan. 

Eigure [TJa) shows a magnetogram of AR 10978 from 


the Michelson Doppler Interferometer (MDI; Scherrer et al. 
1995)1 on board the Solar and Heliospheric Observatory 


(SOHO) on 2007 December 13. The region is isolated, being 
the only numbered active region on the disk at the time. The 
inset in Eigure [^b) shows the soft X-ray emission from the 
core of the region, as observed by Hinode/XRT. The core of 
the active region, where the strongest vertical currents (of crit¬ 
ical importance to NLFEF modeling) are often located, lies 
within the Hinode/SOT-SP scan area, which is demarcated by 
the white boxes in Figures a) and (b). Figures [TJc) and (d) 
illustrate the continuum intensity and the longitudinal mag¬ 
netic field from Hinode/SOT-SP. The region is sufficiently 
small to allow the Hinode held of view to encompass most of 
the magnetic hux. 


In the middle of its disk passage, AR 10978 possessed a 
Hale classihcation of ^ 7 , indicating a high degree of hare pro¬ 
ductivity, and indeed almost 30 soft X-ray Hares are attributed 
to the region in the National Geophysical Data Center GOES 
soft X-ray event list^over the course of its observed history 
on the solar disk. All were relatively small and short events. 
The largest hare was a very short and impulsive C4.5 Hare 
which occurred at 09:39 UT on 2007 December 13, a few 
hours before the Hinode/SOT-SP normal-map scan used here 
was recorded. Additionally, the region produced a Cl.O Hare 


^To download the Level-1 data for this Hinode/SOT-SP scan, please 
see the corresponding entry in the Heliophysics Coverage Registry 
and links therein: http: //www. Imsal. com/hek/hcr ?cmd= 

view-event&event-id=ivo%3A%2F%2Fsot.Imsal.com% 
2FVOEvent%23VOEvent_ObsSP2007-12-13T12%3A18%3A05. 
117.xml 

® Available from ftp://ftp.ngdc.noaa.gov/STP/ 

space-weather/solar-data/solar-features/ 
solar-flares/x-rays/goes/ 
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Fig. 1.— Images of NOAA AR 10978 on 2007 December 13. Panel (a) shows the SOHO/MDI full-disk magnetogram at 12:46 UT, 
obtained within the interval of the Hinode normal-map scan used in this study. The image saturates at ±1000 Mx cm“^. Panel (b) 
shows a logarithmically scaled Hinode/XRT images (Ti/Poly filter) averaged over the scan interval, for context. Representative 
Hinode/SOT-SP data are shown in the two smaller panels (both 162" x 162") at right: panel (c) is the continuum intensity, and 
panel (d) shows the longitudinal magnetic field derived from the Hinode polarization spectra (scaled to ±1500 Mx cm“^). The white 
boxes in panels (a) and (b) correspond to the region represented by panels (c) and (d). The black dashed outline in panel (d) indicates 
the region subsequently remapped into helioplanar coordinates for use in this study. 


within the hour following the completion of the scan. 

In this study, we construct a set of vector magnetograms 
from the chosen Hinode/SOT-SP scan of AR 10978 in order 
to investigate the effects of spatial resolution on the subse¬ 
quent NLFFF extrapolations. Vector magnetograms are pre¬ 
pared at near-Hinode normal-map resolution, and also at spa¬ 
tial resolutions lowered by factors ranging from 2 to 16. The 
northernmost 25% of the Hinode field of view (mostly con¬ 
taining quiet sun) is excluded in order to make the subsequent 
NLFFF extrapolations more computationally tractable. The 
region considered for modeling is contained within the dashed 
box in Figure [^d). 

Because the codes used for NLFFF modeling in this pa¬ 
per assume a Cartesian geometry, it is necessary to project 
the Hinode data onto a regular helioplanar grid. The field of 
view of the Hinode observations is small relative to the solar 
radius, making the effects of curvature small, and as a result 
this remapping is not expected to significantly affect the mod¬ 
eling results. During the remapping process, B is reprojected 
so that it appears as if the active region were located at disk 


center. The resulting set of remapped vector magnetogram^ 
span a region approximately 168 Mmx 124 Mm in size. The 
local helioplanar coordinates are denoted {x, y, z), with x de¬ 
noting solar west direction, y solar north, and 0 the vertical 
direction. 


The procedure followed in preparing the vector magne¬ 
togram data at each resolution level is as follows: (1) the 
Level-1 Hinode/SOT-SP polarization spectra are rebinned by 
the specified factor; (2) a spectral inversion using the High 


Altitude Observatory Milne-Eddington inversion code (Sku- 


manich & Lites| 19871 Leka & Barnes|20T2 1 is performed; (3) 


the 180° disambiguity is resolved using the North West Re¬ 
search Associates “MEO” minimum energy algorithm (Leka 


& Barne^|2012 1 ; (4) helioplanar components of B, and val¬ 
ues for the vertical component of the current density and 
force-free parameter a = yoJz/Bz, are calculated; and (5) 
the resulting values are remapped (interpolated) onto a regu¬ 
lar and uniformly spaced helioplanar grid. 


^The remapped vector magnetogram data, with uncertainties, are available for 
download from http://dx.doi.org/10.7 910 /DVN/KOUAOU 
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Throughout this article, we refer to the various vector mag¬ 
netograms by the factor by which the observed polarization 
spectra are rebinned in step (1). For example, “bin 2” data are 
rebinned by a factor of two in each dimension, “bin 3” data 
are rebinned by a factor of three in each dimension, etc. We 
also retain a “bin 1” dataset that is prepared without any re¬ 
binning. The grid and pixel sizes for these data as a function 
of the bin level are summarized in the first three columns of 
Table □ 

Uncertainties are provided for the components of B, for 
Jz, and for the force-free parameter a, based on propagation 
of uncertainties from the inversion and disambiguation steps. 
The disambiguation uncertainties comprise azimuthal errors 
of 180° assigned to points where multiple trials of the dis¬ 
ambiguation code produced different results with a frequency 
greater than 10% when the optimization procedure used is 
seeded with different sets of random numbers. 

Figure illustrates the reprojected vector magnetogram 
data. The figure shows maps of Bz and Jz for the bin 1 and 
bin 8 data (with the Jz values shown only for points where 
the signal-to-noise ratio in a is greater than 0.25). We have 
set B = 0 at points in the reprojected coordinate system that 
lie outside of the Hinode field of view (such points are located 
at the eastern and western edges of the remapped data). The 
dashed boxes in the figure correspond to the region used for 
comparing NLFFF extrapolations in Section]^ and is slightly 
smaller than the remapped field of view. As expected. Fig¬ 
ure shows that large-scale structures present in Jz at the 
photosphere are retained in the reduced-resolution data (al¬ 
beit with a smaller magnitude), but structures on scales below 
the reduced resolution limit are lost. 

The structure of the magnetic field, which is inferred from 
the rebinned polarization spectra, is affected by the spatial 
resolution, sometimes in non-intuitive ways. Some of the 
trends depend on the underlying structure. To illustrate this 
effect, we use continuum intensity maps from Hinode/SOT- 
SP to segment the boundary data into umbral and penumbral 
regions. Figurej^a) shows a bin 1 continuum intensity image 
with contours outlining umbral and penumbral regions over¬ 
laid. When downsampling the contours from the bin 1 data 
down to the resolution of the bin 16 data and comparing with 
the rebinned intensity image, as shown in Figure [^b), differ¬ 
ences between the downsampled contours and the locations of 
umbral and penumbral pixels are evident. This effect is espe¬ 
cially noticeable in regions in the bin 1 data where small-scale 
features are present, leading to the classification of some pix¬ 
els determined to be in the umbral region in the bin 1 image 
being identified as penumbral in the bin 16 image. 

Figure shows how the magnetogram-averaged absolute 
vertical flux density \Bz\, absolute horizontal flux density 
Bh = {B'l + ByY/'^, \Jz\, and a vary as the spatial reso¬ 
lution changes. In the set of plots, the inferred values of these 
quantities are plotted as a function of spatial resolution in two 
ways, depending on whether the contours used in the segmen¬ 
tation are determined from the rebinned continuum images 
or are determined from downsampling those from the high- 
resolution bin 1 image. The differences between the two pairs 


of curves in {\Bz\) and (Bh) (Figs.|^a) and (b)) as a function 
of resolution are often due to the spatial resolution affecting 
whether pixels are classified as being within the umbral or 
penumbral contours. 

More strikingly, the average vertical current densities 
(I JjI) are seen to decrease as the spatial resolution decreases, 
indicating that the vertical currents often have structure on 
small scales. This effect also affects the trend in a, with both 
umbral and penumbral values getting closer to zero (i.e., po¬ 
tential) as the spatial resolution decreases. The penumbral 
and umbral values of a possess different signs for most bin 
levels, though the mean a value in penumbral areas changes 
sign as the boundary data become less resolved. 

Table lists values of the net magnetic flux $, and the 
components F) (with i = x, y, z) of the magnetic force on 
the coronal volume, as a function of the bin level. The flux is 
calculated by integrating Bz over the magnetogram area, and 
the force components are obtained by performing the integrals 
in |Molodenskii| ( |1969[ l, which are derived from moments of 
Equation Q. The net flux fl) is expressed as a fraction of the 
unsigned magnetic flux $o (the integral of | over the mag¬ 
netogram area), and the force components are in units of Fq, 
the integral of the magnetic pressure over the magnetogram 
area. The net flux is zero if the observed area is in flux bal¬ 
ance, and the force components vanish if the boundary data 
are consistent with the NLFFF model. The values in the table 
indicate that the region is close to flux-balanced. The region 
has relatively small horizontal force components and larger 
vertical forces (cf. |Metcalf et al.|1995| l. 

2.2. Methods and Codes Used 

Many methods of solution of the NLFFF equations exist. 
In the study presented here, implementations of three such 
methods — the optimization, magnetofrictional, and Grad- 
Rubin methods — are applied to the sequence of vector mag¬ 
netogram data described in Section [ZT] Altogether there are 
five codes (one code implementing the optimization method, 
one code implementing the magnetofrictional method, and 
three codes implementing the Grad-Rubin method in differ¬ 
ent ways). 

Each of the five codes solves for a NLEEE on a uniform 
three-dimensional Cartesian grid with equal grid spacing in 
each dimension. The grids used are defined by the helioplanar 
vector magnetogram boundary data at the lower boundary and 
the choice of a vertical height for the grid. The codes differ 
in how they incorporate the vector magnetogram data at the 
lower boundary, and how they use the information provided 
by the uncertainties in the data. The boundary conditions at 
the side and top of the computational domain are not con¬ 
strained by the observational data, and as a consequence the 
different codes handle boundary conditions at these surfaces 
in different ways. In the following subsections, we briefly 
describe each of the specific codes used in this study. 
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Table 1 

Properties'^ of AR 10978 Remapped Vector Magnetogram Data 


Bin Level 

Size (pixels) 

Pixel Scale (Mm) 

$/4>o 

F,/Fo 

Fy/Fo 

F./Fo 

Fo [N] 

1 

1129x837 

0.106 

-0.030 

8.3x10^ 

-1.1x10-2 

-031 

2.4x10'® 

2 

564x418 

0.212 

-0.026 

9.7x10^ 

-1.1x10-2 

-0.37 

2.4x10'® 

3 

375x278 

0.318 

-0.024 

-3.6x10^ 

-1.1x10-2 

-0.36 

2.5x10'® 

4 

282x209 

0.424 

-0.026 

-9.3x10^ 

-9.5x10-2 

-0.36 

2.4x10'® 

6 

187x138 

0.635 

-0.024 

-4.8x10-^ 

-8.2x10-2 

-0.36 

2.5x10'® 

8 

141x104 

0.847 

-0.026 

-5.1x10-^ 

-5.9x10-2 

-0.35 

2.7x10'® 

10 

112x82 

1.06 

-0.031 

-8.8x10-^ 

-5.8x10-2 

-0.35 

2.6x10'® 

12 

93x68 

1.27 

-0.029 

-7.7x10-^ 

-7.3x10-2 

-0.35 

2.8x10'® 

14 

80x58 

1.48 

-0.030 

-1.4x10-2 

-6.0x10-2 

-0.35 

2.8x10'® 

16 

70x52 

1.69 

-0.038 

-1.6x10-2 

-5.0x10-2 

-0.35 

2.9x10'® 


‘^The ratio is the ratio of the net to unsigned flux, and is a measure of the flux imbalance of the magnetogram, and Fx, Fy, 

and Fz ai‘e the magnetogram-integrated integrated Lorentz forces, as determined from the integrals in |Molodenskii|f7969^ . The forces 
are normalized by Fq, the magnetogram-integrated magnetic pressure. 



Fig. 2.— The vertical magnetic field B^, as derived from the Hinode data for AR 10978 and after remapping to helioplanar 
coordinates, at a resolution close to the Hinode observations (bin 1) is shown in panel (a), and the associated vertical current density 
J 2 is shown in panel (b). Panels (c) and (d) illustrate and after rebinning the Hinode data by a factor of eight and remapping 
into helioplanar coordinates. The color bars in each panel indicate the image scale in Mx cm“^ (for B^) and mA (for J^). The 
dashed boxes correspond to the base of the analysis volume used in Sectionj^for comparing the resulting NLFFF solutions. 
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Fig. 3.— Images of continuum intensity for the (a) bin 1 and (b) bin 16 boundary data. The contours in panel (a) indicate the umbral 
(red) and penumbral (purple) boundaries derived from the bin 1 data, while those in panel (b) are derived from the bin 1 data and 
subsequently downsampled to the coarser resolution of the bin 16 data. In panel (b), it is evident that the downsampled umbral and 
penumbral contours do not exactly correspond to the analogous contours that would be drawn using the bin 16 continuum image. 
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Fig. 4.— Variations of the mean values of (a) unsigned vertical magnetic flux density (b) magnitude of the horizontal flux 
density Bh = (B^^ + ByY^^, (c) magnitude of the vertical current density \ Jz\, and (d) force-free parameter a that characterize 
the vector magnetogram data, as averaged over the full area and plotted as a function of spatial resolution. For clarity, we note that 
spatial resolution increases leftward (toward smaller pixel sizes) in each plot. In each panel, the averages over all pixels are shown 
in red, and averages over pixels inside the umbral and penumbral contours are shown in green and blue, respectively. Dashed curves 
indicate that the umbral and penumbral regions were determined from the rebinned data, and solid curves indicate that the regions 
were based on contours downsampled from the bin 1 data. (An example of contours that have been downsampled from bin 1 to 
bin 16 is shown in Fig.j^b).) Error bars indicate standard deviations, except for panel (d) where standard errors are shown. 


















































2.2.1. Optimization Method 


The optimization method is a relaxation scheme that seeks 
to minimize a volume integral such that, if the integral be¬ 
comes zero, the field is divergence- and force-free. In its orig¬ 


inal form (Wheatland et al. 20001, the method proceeds as 
follows. An initial magnetic field is chosen in the computa¬ 
tional volume, and lower boundary values are replaced by the 
required vector boundary conditions B. The field is evolved 
forward using equations that minimize a functional contain¬ 
ing the magnitudes of the Lorentz force and of the divergence 
of B. The process is halted when the field reaches an approx¬ 
imately steady state. 

The optimization solutions discussed here are based on the 
algorithm described in Wiegelmann & Inhester ( 2010) l and 
Wiegelmann et al. (2012|, which includes several modifica¬ 


tions to the original method. As discussed in [Wiegelmann] 
jet al.| ( |200^ , the optimization solutions are found to improve 
when a preprocessing scheme is applied to the vector bound¬ 
ary data. Preprocessing reduces the magnitude of the inte¬ 
grals representing the boundary-integrated magnetic force and 
torque (which necessarily vanish for a NLFFF), subject to 
penalty functions that simultaneously aim to preserve agree¬ 
ment with the observed vector magnetogram data. Spatial 
smoothing is also applied to the boundary data. The spe¬ 
cific preprocessing scheme used here involves four weight¬ 
ing parameters that determine which constraints are most 
closely met. The values used here are /ri = /i 2 = 1 and 
1^3 = M 4 = 10“^, where the weighting parameters are as de¬ 
scribed in Equation (6) in Wiegelmann et al. ( [2006 | l. 

Additionally, modifications of the lower boundary values 
are permitted during the optimization process, implemented 
by including an additional optimization term in the functional 
that takes into account the measurement uncertainties in the 
vector magnetogram data. This term allows more substantial 
changes to the components of the field at points where the 
associated uncertainties are larger. These changes are con¬ 
trolled by a two-dimensional weighting matrix IV (x, y) in an 
optimization integral over the lower boundary, as described in 
Wiegelmann & Inhester (2010|l. 


The modified optimization method also includes three- 
dimensional weighting functions Wf{x, y, z) and Wd{x, y, z) 
in the functionals for the volume-integrated Lorentz force and 
divergence (see Equation (4) of [Wiegelmann et al.||2012| . 
These weights wf and Wd are set to unity in the entire model 
volume, except for finite boundary layers adjacent to the lat¬ 
eral and top boundaries of the computational volume, where 
they smoothly approach zero. This causes the solution ob¬ 
tained by the method to remain fixed at the boundary values 
prescribed by the initial field, and introduces buffer regions at 
the side and top boundaries where the solution field may de¬ 
part from a force-free and divergence-free state. These buffer 
regions are excluded from the analysis in Section]^ 

Another feature of the optimization code used here is a 
grid-refinement scheme. The scheme entails applying the op¬ 
timization code on coarser grids, the solutions of which are 
then used to initialize the optimization algorithm on succes¬ 


sively more refined grids. Typically several refinement lev¬ 
els are used. The series of increasingly more refined grids is 
started by using a potential field to initialize the coarsest grid. 
This scheme has been shown in earlier studies to improve the 
quality of the resulting solutions to the NLEEE model, and 
decreases the running time of the full calculation. 


2.2.2. Magnetofrictional Method 


Magnetofrictional codes evolve the magnetic induction 
equation using a velocity field that advances the solution to 
a more force-free state (e.g., Chodura & Schlueter[1981|l. The 


magnetofrictional code used here is described in Valori et al. 


( [2007 2010| l, and its application to solar data is presented in 
[Valori et~ ( [2012b| l. The magnetofrictional method uses the 
full vector field over the entire lower boundary as boundary 
condition, taking as input preprocessed vector magnetogram 
data. Both the horizontal and vertical components of B may 
be adjusted during the preprocessing stage to make the bound¬ 
ary data more compatible with the force-free model, although 
changes in the boundary data are constra ined to be within cer¬ 
tain limits (Fuhrmann et al. 2007 2011 


As the field is evolved forward in time, the boundary con¬ 
ditions at the side and top boundaries are constructed by 
the code at each iteration of the magnetofrictional relaxation 
method using the neighboring volume values of the field at 
the previous iteration ( [Valori et al.[|2007'j l. The construction 
uses first order polynomial interpolation, and requires that the 
solenoidal condition is met at the boundary points, and that 
the Lorentz force falls to zero at ghost points just exterior to 
the volume. 


A grid-refinement strategy analogous to that used in the 
optimization code is also used here. The domain having 
the coarsest refinement is initialized using a potential mag¬ 
netic field with values of Bz matching the initial prepro¬ 
cessed lower boundary values, and subsequent, more refined 
domains are initialized using the resulting fields from the less 
refined solutions. 


Eor this particular active region, the magnetofrictional ex¬ 
trapolations of the full field of view at all tested resolutions 
developed abnormally strong but stable large-scale currents 
near the northern lateral boundary. Such behavior is unusual, 
and remained even after tuning the parameters of the algo¬ 
rithm. As a result, additional cropping of the magnetogram 
at the northern side was performed before running the mag¬ 
netofrictional models shown here. In all but the case using 
bin 3 boundary data, the cropped magnetograms allowed the 
preprocessing algorithm to produce more compatible (more 
force-free) boundary conditions for the extrapolation code, 
and the resulting extrapolated fields did not develop any un¬ 
usual current at the northern lateral boundary. Eor the bin 3 
case, even this additional cropping did not yield a consistent 
preprocessed boundary condition, and consequently this case 
is omitted from the analyses presented in Section]^ 


®Both the magnetofrictional and optimization codes apply preprocessing algo¬ 
rithms to the boundary data, however the preprocessing algorithms used with 
these codes are distinct and unrelated. 


8 






















































2.2.3. Grad-Rubin Method 


XTRAPOL code 


Three different implementations of the Grad-Rubin method 
are used. These are the CFIT code ( Wheatland 2007| l; the 
XTRAPOL code (as described in |Amari et al.|200^ and with 
the boundary co nditions given in |Amari & Aly 2010) l; and 
the FEMQ code ( |Amari et al.|2006| . A brief account of the 
specific implementations and boundary conditions is given 
here. 


Grad-Rubin methods take as boundary conditions the nor¬ 
mal component of the field in the boundary and the value of 
the force-free parameter a over one polarity of the field in the 
boundary. Since there are two choices of polarity there are 
two solutions, which are in general different if the boundary 
data are inconsistent with the force-free model. We present 
both solutions for all Grad-Rubin codes in this paper, using P 
and N to denote solutions with a chosen from positive- and 
negative-polarity regions in the lower boundary, respectively. 


The XTRAPOL code solves a mixed elliptic-hyperbolic 
boundary-value problem for a and B, which is mathemat¬ 
ically well posed ( [Boulmezaoud & Amari 2000) l. It uses a 
finite-difference approach with a representation of B based 
on a vector potential A (defined with a convenient gauge) on 
a staggered mesh. This ensures that the divergence operator 
remains exactly on the kernel of the curl operator (and thus 
V ■ B = 0 to rounding errors), independently of the mesh 
resolution. It solves iteratively the elliptic problem through 
a positive definite linear system, and the hyperbolic problem 
by transporting a along the characteristics (the magnetic field 
lines), imposing a originating from only one or the other 
polarity. The code uses the MPI library. The solution can 
be provided for both balanced or non-balanced photospheric 
magnetic flux, in which case field lines can intersect other 
external boundaries. SeejAmari et al.|(|2006|l for more details. 


CFIT code 


FEMQ code 


Grad-Rubin methods involve updates to both B and J at each 
iteration, involving the solution of linear hyperbolic and ellip¬ 
tic partial differential equations, respectively. The CFIT code 
performs the update to J by propagating lower boundary val¬ 
ues of a along numerically traced field lines. The update to B 
is achieved by solving the Poisson equation for the vector po¬ 
tential for the field using a two-dimensional Fourier Trans¬ 
form method, with transforms in the x and y directions. Solu¬ 
tions from the code are correspondingly periodic in x and y. 

With the CFIT code, the lower boundary conditions on the 
field are the values of Bz from the vector magnetogram data. 
Boundary conditions on a are chosen as follows. Values of a 
calculated using Equation 0 are used at a given boundary 
point provided the signal-to-noise ratio is at least 5%, based 
on the uncertainties provided with the data, and provided the 
boundary point does not exhibit a large localized spike in a. 
Boundary points not meeting these criteria are assigned values 
a = 0. This procedure is referred to as “censoring” of the 
boundary data. 

Additionally, during the Grad-Rubin iteration procedure, 
points in the domain that are threaded by field lines that inter¬ 
sect the top or side boundaries are assigned a = 0 to ensure 
V • J = 0 globally. As a result, boundary points in such open- 
field regions are treated as having a = 0, however the regions 
of open and closed field lines generally change as the calcula¬ 
tion proceeds, and consequently boundary points that are thus 
censored at one iteration may not be censored at another. 


The iteration is initiated with a potential field in the vol¬ 
ume with values of Bz matching the vector magnetogram at 
the lower boundary values, obtained by a two-dimensional 
Eourier transform solution. The CEIT code is parallelized us¬ 
ing both the OpenMP ( Chandra et al.|200T) l and MPI (Gropp 


|et al.|I999|) standards. Eurther details on the solution method 


and the handling of boundary conditions are given in Wheat- 
[Iandl ( |2007l l. 


The EEMQ code solves the same well posed boundary-value 
problem as XTRAPOL. However EEMQ uses a finite-element 
approach, and works directly with B using a least-squares 
approach that minimizes the divergence of B. The hyperbolic 
equations at each iteration can be solved either using the same 
method as in XTRAPOL or by solving a non-positive-definite 
linear system for a by an iterative method. More details can 
be found in [Amari et al.| ( |2006| . 


3. Results 


This section presents quantitative comparisons amongst 
the extrapolated solution fields and characterizes the effects 
of the spatial resolution of the boundary data on the solutions. 
We analyze the energy and helicity of the solutions and dis¬ 
cuss the departure of the lower boundary field values in the 
resultant solutions from the vector magnetogram values. The 
magnetic energy analysis includes a Helmholtz decomposi¬ 
tion of the solution fields into solenoidal and non-solenoidal 
components, following the procedure of Valori et al. (2013 i, 
and a discussion of the effects of the non-solenoidal compo¬ 
nent on the resulting energies. Because physical magnetic 
fields are solenoidal (due to V • B = 0), comparing the non- 
solenoidal and solenoidal components provides a check on the 
consistency of the solutions. 


Performing a NLEEP extrapolation using the bin 1 data 
remains a challenge for the codes and for computer hard¬ 
ware. An extrapolated field B corresponding to these data 
and represented by a single-precision, floating-point, three- 
dimensional array possesses of order «1000^ points and 
requires «12 GB of computer memory. The codes as written 
require multiple copies of this and similar three-dimensional 
arrays to be stored simultaneously in computer memory. As 
a result, memory use is an issue. A more constraining factor 
is the scaling of the time to completion, taken to be at best 
cx for serial calculation on a grid with points (e.g., 
|Wheatland|2007|). A number of the methods employ parallel 
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implementations to improve performance, but the total com¬ 
putational time remains a problem for large grids. As a result 
of these constraints, we do not present the results obtained for 
the bin 1 data. 

Results are presented for bin levels 2-16 for all five codes 
(the optimization and magnetofrictional method codes, and 
the CFIT, XTRAPOL and FEMQ Grad-Rubin codes), with 
the exception of the magnetofrictional solution using the bin 3 
boundary data. For the three Grad-Rubin codes, results are 
presented for both the P and N cases. Altogether, there are 
71 solution data cubes 0 

For a given resolution level, the five codes employ domains 
that are slightly different. To standardize the comparisons 
across codes and resolutions, we will use a fixed analysis vol¬ 
ume V, chosen to be the largest common physical volume for 
all solutions across codes and resolutions. This analysis vol¬ 
ume is 136x77x107 Mm^ and covers approximately the re¬ 
gion from X = 14.3 Mm to a: = 150.0 Mm, from y = 14.3 Mm 
to y = 92.8 Mm, and from z = 0 Mm to z =107 Mm, where 
the origin is located at the lower left-hand corner of the input 
magnetograms0 The footprint of the analysis volume V on 
the photosphere is indicated by the dashed lines in the vector 
magnetogram data shown in Figure]^ 

Figure [^presents visualizations of field lines within V for 
the solutions for the different methods for the finest resolu¬ 
tion level for which there are solution data. Qualitatively, it 
is evident that there are some differences between the results 
obtained with the different methods or codes at this resolu¬ 
tion. Such variation across the solution methods is in com¬ 
mon with our earlier studies (e.g., Schrijver^etalJ2006yMet 


calf et al!]|2008[ Schrijver et al.||^08[ De Rosa et al.||2009 1. 
When examining the differences in solution data across the 
range of resolution levels, as calculated by any of the codes, 
the initial qualitative impression is that there are only minor 
differences in (for example) the detailed field-line trajecto¬ 
ries. The shape of the field-line bundles originating in differ¬ 
ent locations across the lower boundary of the analysis region 
appears similar for all resolutions, and the boundaries sepa¬ 
rating field lines that leave the analysis volume versus those 
that are contained within the domain appear similar. In the 
sections that follow, we show that more significant variations 
exist amongst the solutions, both across resolution levels and 
across methods and codes, when quantitative comparisons are 
performed. 


3.1. Magnetic field energy, helicity and extrapolation 
metrics 

In this section, the magnetic energies E, the relative mag¬ 
netic helicities Hm of the extrapolated magnetic fields, and 
additional extrapolation metrics are discussed. The variations 
amongst the energies and helicities of the extrapolated fields 
depend on many factors, including the effects of resolution, 
the boundary conditions taken as input, and the handling of 
boundary conditions by the various codes. An additional con¬ 
sideration is that, for each resolution level, the same numeri¬ 
cal grid is used for the calculation of these quantities, regard¬ 
less of code and method, with interpolation onto this grid per¬ 
formed when necessary (such as when analyzing results from 
XTRAPOL, in which variables are offset with each other on 
a staggered mesh). Centered, second-order finite differences 
are used to calculate the needed derivatives, and this differ¬ 
encing scheme may also be different than the one used by the 
various codes (e.g., the magnetofrictional code uses fourth- 
order differences). 

3.1.1. Magnetic Energy 

Figure illustrates the magnetic energy E, the magnetic 
energy of the corresponding reference potential field Eq, the 
free energy Ef = E — Eq, and the ratio E/Eq = 1-1- Ef /Eq 
as a function of spatial resolution for each of the extrapo¬ 
lations within the common analysis volume V. The refer¬ 
ence potential fields are computed numerically by solving 
the Laplace equation for the potential with Neumann bound¬ 
ary conditions based on the normal component of B on 
all six boundaries of the analysis volume (which we denote 


®The solution volumes used for analysis are available for download from 
http://dx.doi.org/10.7 910/DVN/7ZGD9P 


Bn\dv)- Following Thomson’s theorem (e.g., Jackson 1999 i, 
the potential field uniquely represents the field with the mini¬ 
mum energy for any divergence-free field where Bn matches 
Bn\dv- As a consequence of each solution field having its 
own Bn\av^ there are separate reference potential fields for 
each of the solution fields. 

The energy metrics are listed in Table|^ Taken as a group, 
the energies of the solution fields range from 1.08 x 1026 J to 
1.50 X1026 J, with a mean value of 1.25 x lO^® J and a standard 
deviation of 10%. In terms of free energy Ef, the solutions 
range from near-potential fields having free energy 2% above 
Eq (both P and N models of XTRAPOL using the bin 16 
boundary data) to a case where the free energy is 24% of the 
associated potential field (optimization using the bin 2 bound¬ 
ary data). In physical units, the free energy is in the range 
0.3-2.9x1026 J. 


^The indices of the analysis sub-volume V for each resolution, relative to 
the boundary data, are as follows: bin 2: [48:504, 48:312, 0:360] (size 
457x265x361); bin 3: [32:336, 32:208, 0:240] (size 305x 177x241); bin 4: 
[24:252, 24:156,0:180] (size 229x 133x 181); bin 6: [16:167,16:103,0:119] 
(size 152x88x120); bin 8: [12:126, 12:78, 0:90] (size 115x67x91); 
bin 10: [10:100, 10:62, 0:71] (size 91x53x72); bin 12: [8:84, 8:52, 0:60] 
(size 77x45x61); bin 14: [7:71, 7:44, 0:51] (size 65x38x52); bin 16: 
[6:63, 6:39, 0:45] (size 58x34x46). 


Figures |^c) and (d) show that the solutions obtained with 
the magnetofrictional and optimization methods have higher 
free energies than those calculated using Grad-Rubin methods 
(CFIT, XTRAPOL, and FEMQ), both in absolute (E - Eq) as 
well as in relative {E/Eq) terms. When more highly resolved 
boundary data are used, all methods return solutions that trend 
toward higher energies, with the exception of the magnetofric¬ 
tional method. Eigure illustrates that the free energies for 
the magnetofrictional method range from 10% to 19%, with 
higher resolution cases generally having less free energy. In 
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Table 2 

NLFFF Extrapolation Metrics^ for AR 10978 


Method/Code 

Bin 

E [10^6 J] 

EjEo 

Hm [ 10^6 Wb 2 ] 

(CWsine) 

(l/r|)[xl 0 ^] 


Optimization 

16 

1.41 

1.12 

0.66 

0.13 

19. 

0.38 


14 

1.38 

1.14 

0.03 

0.14 

18. 

0.37 


12 

1.43 

1.18 

-0.16 

0.12 

16. 

0.39 


10 

1.39 

1.15 

-0.24 

0.14 

10 . 

0.33 


8 

1.46 

1.18 

0.47 

0.11 

6.1 

0.23 


6 

1.43 

1.18 

-0.63 

0.11 

3.7 

0.21 


4 

1.48 

1.21 

0.24 

0.10 

2.2 

0.18 


3 

1.47 

1.20 

-0.92 

0.13 

1.7 

0.19 


2 

1.50 

1.24 

0.04 

0.10 

1.1 

0.16 

Magnetofrictional 

16 

1.37 

1.16 

1.29 

0.30 

82. 

0.51 


14 

1.22 

1.15 

3.79 

0.30 

71. 

0.51 


12 

1.23 

1.17 

2.62 

0.32 

52. 

0.53 


10 

1.26 

1.19 

3.65 

0.25 

46. 

0.51 


8 

1.28 

1.14 

4.89 

0.31 

33. 

0.34 


6 

1.11 

1.09 

3.17 

0.29 

18. 

0.51 


4 

1.23 

1.12 

3.96 

0.27 

13. 

0.26 


2 

1.08 

1.10 

1.80 

0.29 

13. 

0.34 

CFIT (P / N) 

16 

1 . 12 / 1.12 

1.04/1.03 

3.55 / 3.43 

0.35 / 0.40 

11 . / 12 . 

0.08 / 0.09 


14 

1.09/1.08 

1.06/1.04 

4.14/3.38 

0.32/0.38 

9.9/11. 

0.07/0.10 


12 

1.11 / 1.11 

1.05/1.04 

3.56/3.95 

0.32/0.36 

7.6/ 9.8 

0.06/0.10 


10 

1.09/1.09 

1.05/1.05 

3.51/4.23 

0.32/0.34 

6.5/ 8.6 

0.08/0.11 


8 

1 . 12 / 1.12 

1.05/1.05 

3.32/4.58 

0.32/0.31 

6.4/ 7.3 

0.12/0.13 


6 

1.11 / 1.09 

1.07/1.05 

4.05 / 3.68 

0.27 / 0.29 

4.7/ 6.1 

0.10/0.14 


4 

1.19/1.16 

1.11 / 1.10 

5.21/3.88 

0.28/0.30 

5.0/ 6.1 

0.21/0.24 


3 

1.11 / 1.11 

1.06/1.06 

2.45 / 2.95 

0.27 / 0.26 

4.0/ 4.8 

0.21/0.24 


2 

1 . 12 / 1.10 

1.06/1.05 

4.47 / 2.32 

0.24 / 0.25 

2.1/ 3.4 

0.19/0.25 

XTRAPOL (P / N) 

16 

1.29/1.29 

1 . 02 / 1.02 

2.68/2.30 

0.26 / 0.25 

7.9 /7.8 

0.03 / 0.03 


14 

1.25/1.25 

1.03/1.03 

3.01/3.01 

0.25 / 0.24 

6.7 /6.5 

0.03 / 0.02 


12 

1.28/1.27 

1.03/1.03 

3.21/2.99 

0.24 / 0.22 

6.0 /5.7 

0.03 / 0.02 


10 

1.27/1.26 

1.04/1.03 

3.99/2.75 

0.21/0.19 

5.1 /4.7 

0.03 / 0.02 


8 

1.29/1.28 

1.04/1.04 

3.47 / 3.24 

0.21/0.18 

3.1 /3.0 

0.02 / 0.02 


6 

1.27/1.26 

1.05/1.04 

3.45 / 2.72 

0.17/0.16 

2.2 / 2.1 

0.02 / 0.02 


4 

1.29/1.28 

1.05/1.05 

3.52/4.13 

0.13/0.16 

1.2 / 1.2 

0.02 / 0.04 


3 

1.28/1.27 

1.05/1.04 

3.40/2.83 

0 . 12 / 0.12 

0.83/0.77 

0.02 / 0.02 


2 

1.28/1.28 

1.05/1.05 

3.46 / 2.24 

0.14/0.13 

0.46 / 0.40 

0.08 / 0.06 

FEMQ (P / N) 

16 

1.29/1.30 

1.04/1.03 

2.87/2.50 

0.28 / 0.25 

6.4 /6.9 

0.02 / 0.02 


14 

1.25/1.25 

1.03/1.03 

3.32/3.14 

0.26 / 0.26 

5.2 /5.3 

0.02 / 0.02 


12 

1.28/1.28 

1.03/1.03 

3.24 / 3.06 

0.25 / 0.24 

4.5 /4.6 

0.02 / 0.02 


10 

1.27/1.26 

1.04/1.04 

3.87/2.70 

0.22 / 0.20 

3.5 /3.5 

0.02 / 0.02 


8 

1.29/1.28 

1.05/1.04 

3.93/3.50 

0.21/0.19 

2.8 / 2.8 

0.02 / 0.02 


6 

1.27/1.26 

1.05/1.04 

3.65/2.65 

0.17/0.15 

1.8 / 1.7 

0.02 / 0.01 


4 

1.30/1.28 

1.05/1.05 

3.66/3.87 

0.12/0.15 

0.99 / 0.98 

0.01/0.04 


3 

1.29/1.28 

1.05/1.04 

3.34/2.79 

0 . 11 / 0.11 

0.69 / 0.66 

0.04 / 0.03 


2 

1.29/1.28 

1.06/1.05 

3.87/2.97 

0 . 10 / 0.11 

0.42 / 0.38 

0.07 / 0.07 


“All metrics are defined in Section lrT] 


11 











(a) CFIT P 


(d) XTRAPOL P 


(g) FEMQ P 



(e) XTRAPOL N 


(h) FEMQ N 



(c) optimization 


(f) magneto- 
frictional 


(i) potential 


Fig. 5.— Field lines for all solutions using the bin 2 boundary data. The panels show results using the (a) CFIT Grad-Rubin 
(P solution), (b) CFIT Grad-Rubin (N solution), (c) optimization, (d) XTRAPOL Grad-Rubin {P solution), (e) XTRAPOL Grad- 
Rubin (N solution), (f) magnetofrictional, (g) FEMQ Grad-Rubin (P solution), and (h) FEMQ Grad-Rubin (N solution) codes. Eor 
comparison, field lines from a potential field matching the values of provided to the modelers are shown in panel (i). Eield lines 
are plotted in black, green, and magenta: black field lines close within the volume, whilst the magenta and green field lines are open, 
and originate in different polarities. The starting points for the field line trajectories form a regularly-spaced grid at z = 0. 


contrast, the free energies from the optimization method in¬ 
crease as the resolution increases from 12% to 24%. The 
free energies in the Grad-Rubin codes are also seen to double 
over the range of resolutions, becoming greater for the more 
spatially resolved cases (as with the optimization method), 
though this variation corresponds to changes of a few percent 
of the total magnetic energy. Within the set of Grad-Rubin 
solutions, the relative free energy values are slightly higher at 
all resolutions for CEIT than for EEMQ and XTRAPOL. Tree 
energies for the P solutions are equal or marginally higher 
than those from the N solutions. At each resolution value, 
the spread in free energies for the Grad-Rubin codes appears 
small, with the exception of the CEIT calculations using the 
bin 4 boundary data (addressed further in the next subsection). 

Two aspects of the calculation methods may explain the 


spreads in total and free energies. Eirst, the larger spread 
in energies from the optimization and magnetofrictional so¬ 
lutions may be attributed in part to their divergences, as eval¬ 
uated numerically, being greater than the Grad-Rubin meth¬ 
ods (which is expected on the basis of past results, see, e.g., 
Schrijver et al.|20()6| 2008[ De Rosa et al. 2009| l. These meth¬ 


ods introduce a departure from V • i? = 0 upon initialization, 
and then seek to minimize this error during iteration. In con¬ 
trast, the Grad-Rubin implementations employ schemes that 
are initialized with divergence-free fields, and aim to provide 
divergence-free fields at all iterations. (We further address the 
contribution from residual non-solenoidal field components to 
the total energy in the next subsection.) 

Second, the higher energies of the magnetofrictional and 
optimization solutions may be due in part to the preprocess- 
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Fig. 6. — (a) Total magnetic energy E, (b) energy of the associated potential field Eq, (c) free energy Ef = E — Eq, and (d) ratio 
E/Eq, as a function of spatial resolution for the set of NLFFF calculations. Please note that resolution increases to the left in each 
plot (i.e., smaller pixel sizes and smaller bin factors correspond to more highly resolved boundary data). 
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Fig. 7.— The decomposition of the total energy into solenoidal (i7o,s and iJj s) and non-solenoidal (ii^o.ns, £’j,ns, and |£^mix|) con¬ 
tributions, for the different solutions. The individual panels show results for the (a) magnetofrictional, (b) optimization, (c) CFIT, 
(d) XTRAPOL, and (e) FEMQ codes. The different contributions are shown as a function of resolution. The tilde over each term 
indicates normalization to the corresponding total energy E. Table in Appendix [A| lists the numerical values of all contributions 
for each solution. 


ing used with these methods. As was the case in Metcalf et al. 
( 2008 | l and [Schrijver et al. (2008 1 , preprocessing the bound¬ 
ary data to be more compatible with the force-free assumption 
can result in solution fields having greater free energies. All 
Grad-Rubin implementations presented here remove spikes in 
a (many of which likely result from poorly constrained values 
of Bx and By), but do not otherwise preprocess the bound¬ 
ary data. The degree to which the codes necessarily alter the 


boundary data are discussed further in Section 3.2 


3.1.2. Energy Decomposition 

Accurate numerical solutions to the force-free model 
should be solenoidal, to within numerical errors. However, 
the methods of solution to the equations, and the inconsis¬ 
tency of the boundary data with the model, can lead to sig¬ 
nificant departures from solenoidality. |Valori et al.| ( |2013| l 
provide a method for quantifying how much a non-zero diver¬ 
gence in a NLFFF solution affects the accuracy of the estimate 


for the magnetic field energy. The Valori et al. (20131 method 


is an application of Thomson’s theorem. The magnetic field 


is decomposed into a potential and a current carrying part, 
and each of these is split into solenoidal and non-solenoidal 
components via Helmholtz decomposition. The energies of 
the components are then compared. In this section, we sum¬ 
marize the decomposition as applied to each of the NLFFF 
calculations, leaving the details to Appendix [A| 

The magnetic energy within V may be written 

E = — [ B^dV 

2/ro Jv 

= f?0,s + + f^'O.ns + Es „s + F/mix , (6) 


where i?o,s and Sj s are the energies of the potential and 
current-carrying solenoidal components, i?o,ns and E'j ns are 
those of the non-solenoidal components, and E^nix is a non- 


solenoidal mixed term (see Equations (7) and (8) in Valori 


et al. 2013 for the expressions for the energies in terms of 


the field components). The energies in Equation (j^ are posi¬ 
tive by definition, with the exception of E^^ix- For a perfectly 
divergence-free field, i7o,s = Eq and E]^s = Ef = E — Eq, 
while F^o.ns = F'j.ns = F'mix = 0. In practice, a NLEEE calcu¬ 
lated numerically will not be perfectly solenoidal, and So,ns, 
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and Smix will instead be finite. Comparing the non- 
solenoidal components with the free energy can be used to 
gauge the reliability and uncertainty of the free energy deter¬ 
mination. 

The values of each term in Equation (j^ for each method 
and spatial resolution are plotted in Figure For the Grad- 
Rubin methods, the results are plotted for both the N and 
P solutions using identical (and often overlapping) sym¬ 
bolism. Figure shows that the energies contained in the 
non-solenoidal components are generally lower than the 
solenoidal portions. The primary exceptions to this trend are 
with the optimization solutions, where the magnitudes of the 
mixed terms jF'mixl for all resolution levels are of the same or¬ 
der as the free energy i?j j. The energy of the non-solenoidal 
component of the potential field Eq ^s is smaller by at least a 
factor of 10^ in all cases, and in almost all cases decreases for 
higher resolutions. i?o,ns is generally smaller than E^^s and 
Ifi'mixl, which indicates that the potential component of each 
field is generally closer to a solenoidal state than the current- 
carrying component of the field. Of Ej „s and |ii^mix|i the latter 
is usually larger. The elevated values for E'j ns and |i?tnix| for 
the CFIT solutions using the bin 4 boundary data contribute 
to the outlying points (i.e., points that appear to not follow the 
apparent trend) in the total energy plots shown Figure]^ 


3.1.3. Relative Magnetic Helicity 


For each calculation, the relative magnetic helicity Elm in 


the analysis volume V is computed using the method in Val- 


ori et al.|P012a||. The resulting values for all solutions are 
fablel^a ~ 


listed in Table|^nd plotted in Figure]^ We find a wide range 
of values for Hm, with no clear trend that depends on spatial 
resolution but significant dependence upon the extrapolation 
method employed. Values of Hm across all methods and res¬ 
olutions are predominantly positive, except for a few of the 
optimization method solutions. 


The solutions calculated using the three Grad-Rubin 
implementations possess values of Hm that lie between 
2x10^^ Wb^ and about 5x10^® Wb^. The variations with 
resolution within one implementation are found to be of a 
similar order as the variations between the implementations. 
The XTRAPOL and FEMQ codes return the most consistent 
results between different resolutions, and between P and N 
solutions, with a mean relative helicity (over all XTRAPOL 
and FEMQ solutions) of (3.2±0.5)x 10^® Wbl The CFIT 
code solutions exhibit slightly greater variation with spatial 
resolution; the average helicity over all CFIT solutions is 
(3.7±0.7 )x 102'= Wb^. 


The values of Hm from the magnetofrictional solutions 
trend higher and then lower as the resolution changes, with 
a mean of (3.1±1.2)x 10^^ Wb^, in somewhat general agree¬ 
ment with the values of Hm found from the Grad-Rubin so¬ 
lutions. The optimization method solutions have the lowest 
Hm values — including several that are negative — with 
a mean and standard deviation of (—0.1±0.5)x 10^^ Wb^. 
The larger fractional variation in Hm values obtained by the 
magnetofrictional and optimization methods compared with 
the Grad-Rubin codes may be partly explained by the larger 


Binning factor 


234 6 8 10121416 



Fig. 8.— The relative magnetic helicity Hm as a function of 
resolution, for the different solutions. 
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non-solenoidal errors in these solutions, as outlined in Sec¬ 
tion 3.1.2 Additionally, preprocessing has been applied in¬ 


dependently to the magnetograms at different resolutions for 
these methods, which may lead to greater variety in the values 
of Hm for the extrapolated fields. 


3.1.4. Additional Extrapolation Metrics 

The last three columns of Table list domain-averaged 
metrics that have been used in earlier studies. These are 
(CWsin 6 *), the mean sine of the angle 0 between J and 
B at each point, weighted by J; the fractional flux ratio 
(|/i|), where \ fi\ = |(V • B).i\/Ax) and where Ax 
is the grid spacing; and i^, the average of the magnitude of 
the Lorentz force to the sum of the magnitudes of its con¬ 
stituent pressure and tension components (see Equation (10) 
in |Malanushenko et al.|2014| ). 

The (CW sin 6) and ^ metrics are measures of how force- 
free a solution field is, both of which vanish for perfectly 
force-free fields. The Grad-Rubin solutions show similar (al¬ 
beit opposite) trends for both metrics, with (CW sin 9) de¬ 
creasing as more highly resolved boundary data are used and 
^ mostly showing the reverse behavior. No significant dif¬ 
ference is apparent between the P and N solutions. For the 
magnetofrictional method, neither (CW sin 9) nor ^ exhibit 
well defined trends, though there is a tendency for both met¬ 
rics to indicate more force-free solutions as the resolution of 
the boundary data is increased. The solutions calculated by 
the optimization method do show a clear trend for with the 
more highly resolved fields having lower values of however 
any trend for the (CW sin 9) metric is less clear, though the 
general tendency is for more force-free solutions when higher 
resolution boundary data are used. The {\fi\) metric is a mea¬ 
sure of how divergence-free a solution field is, vanishing for 
divergence-free fields. For all methods, this metric decreases 
as more highly resolved boundary data are used. 

3.2. Changes to the Boundary Data 

All of the solution methods introduce significant changes 
to the vector magnetogram boundary conditions, either ini¬ 
tially, and/or during solution of the NLFFF equations. The 
magnitude of the changes are based on the uncertainties asso¬ 
ciated with the remapped vector magnetogram data, and thus 
tend to be more significant in weak-field regions than in re¬ 
gions of stronger field. Here, we characterize these changes 
by directly comparing the values of B on the lower boundary 
of the resultant NLFFF solutions to the values of B supplied 
by the vector magnetogram boundary data. 

For each solution we define the change A in the vector 
field at z = 0 as 

A = B- (7) 

comparing the solution B with the disambiguated vector 
magnetogram value B'^^. The normalized change A ac¬ 
counts for the uncertainties ai in each component i, such that 

~ ~ ~ VM 

A = B-B , (8) 


where Bi = Bijai and = B'f^ j a i, for each compo¬ 
nent i = {x, y, z}. We consider separately the magnitudes of 
the vertical components of the changes (A^ and A^), and the 
magnitudes of the horizontal cornponents^f the changes, de¬ 
fined by Ah = |(Aa;,Ay)| and Ah = \{A^, Ay)\/V2, with 
the factor of v /2 introduced so that the expected value for nor¬ 
malized changes in the horizontal field, taking into account its 
uncertainties, is unity. 

Table l^lists root-mean-square (rms) values for the magni¬ 
tudes of the vertical and horizontal components of the changes 
in the boundary conditions defined by Equations (0 and ([^ 
for all solutions, for boundary points within the footprint of 
the analysis volume V. In calculating the rms values for the 
normalized changes, boundary points with assigned uncer¬ 
tainties less than 1 G have the uncertainty replaced by a value 
of 1 G. This step prevents the small uncertainties biasing the 
rms values. 

The table shows that all methods introduce changes in the 
horizontal field boundary values substantially larger than the 
uncertainties assigned to the vector magnetogram data. In ab¬ 
solute units, the magnitude A™® of the changes in the hor¬ 
izontal field are largely similar (several hundred Gauss) for 
the different methods. The magnetofrictional method intro¬ 
duces the smallest changes, with rms values of order 230 G 
to 240 G, and the other methods introduce changes which are 
larger by about 50%. In a normalized sense, any resolution 
dependence on A™* appears weak. The Grad-Rubin methods 
show no substantial difference in the values of A™'* for the P 
and N solutions. 


In addition, the optimization and magnetofrictional meth¬ 
ods introduce changes in the vertical field, as is evident from 
Table As with the horizontal field changes, these also ex¬ 
ceed the uncertainties. In absolute units. A™® for the op¬ 
timization method decreases from about 160 G when using 
more coarsely resolved boundary data to about 30 G when 
more finely resolved data are used. Similarly, there is a clear 
trend in the normalized change A™% which is found to de¬ 
crease from a factor of 7.89 when the coarsest boundary data 
are used to 1.17 for the data with the highest resolution. The 
magnetofrictional method consistently makes adjustments be¬ 
tween 80 G and 85 G for all resolution levels. As with the op¬ 
timization method, the normalized changes A™* for the mag¬ 
netofrictional method trend smaller as more highly resolved 
data are used, ranging from 4.77 for the bin 16 data to 3.43 
for the bin 2 data. By design, the Grad-Rubin methods (CFIT, 
XTRAPOL, and FEMQ) preserve the values of the vertical 
field on the boundary and thus there are no changes to the 
vertical components. 


For the magnetofrictional method, the changes in the 
boundary values are introduced by the preprocessing applied 
to the boundary data, prior to the NLFFF calculation. The 
changes are comparable to the maximum values allowed by 
the preprocessing (see Section 2.2.2 1 . The specific method of 
preprocessing used constrains the point-wise changes in the 
boundary values to ±100 G (for Bz), and ±150 G (for B^ 
and By). The maximum/minimum values correspond to rms 
values A^ = 100 G and Ah « 280 G. For the optimization 
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Table 3 

Changes to vector magnetogram boundary conditions 


Method/Code 

Bin 

Ar [G] 

Ar [G] 

A rms 

A rms 

Optimization 

16 

161. 

307. 

7.89 

22.2 


14 

117. 

281. 

5.96 

21.3 


12 

134. 

319. 

6.81 

23.8 


10 

103. 

286. 

4.93 

21.6 


8 

117. 

322. 

5.32 

23.7 


6 

117. 

289. 

5.32 

21.3 


4 

80.9 

337. 

3.38 

24.1 


3 

72.4 

301. 

2.88 

20.5 


2 

33.4 

365. 

1.17 

22.1 

Magnetofrictional 

16 

80.0 

225. 

4.77 

14.7 


14 

81.2 

228. 

4.76 

15.1 


12 

80.6 

226. 

4.54 

14.8 


10 

82.2 

229. 

4.47 

14.6 


8 

82.4 

231. 

4.41 

14.1 


6 

82.9 

234. 

4.22 

14.0 


4 

83.9 

237. 

3.96 

13.4 


2 

84.9 

239. 

3.43 

11.7 

CFIT (P/ N) 

16 


348. / 351. 


21.2/21.4 


14 


342. / 349. 


21 .6/22.1 


12 


354. / 361. 


22.3 / 22.8 


10 


370. / 369. 


21.7/22.3 


8 


381. / 380. 


22.2 / 22.8 


6 


384. / 379. 


22.3 / 22.7 


4 


395. / 387. 


22 .2/22.1 


3 


401. /401. 


22.4 / 22.2 


2 


408. / 409. 


20.8 / 20.7 

XTRAPOL (P/N) 

16 


304. / 300. 


21.4/22.1 


14 


306. / 309. 


21.9/22.9 


12 


308. / 315. 


22.3/23.5 


10 


317. /318. 


22.2/23.0 


8 


323. / 334. 


22.3/23.2 


6 


332. / 342. 


22.3/23.5 


4 


349. / 361. 


22.7 / 24.0 


3 


362. / 370. 


22.2/23.4 


2 


382. / 387. 


20.9/21.6 

FEMQ (P/N) 

16 


322. / 325. 


22.0 / 22.6 


14 


327. / 333. 


22.5/23.3 


12 


328. / 339. 


30.4 / 24.0 


10 


335. / 341. 


22.8/23.7 


8 


339. / 354. 


22.7/23.9 


6 


352. / 362. 


22.9 / 24.2 


4 


366. / 378. 


23.2 / 24.5 


3 


375. / 382. 


22.7/23.7 


2 


386. / 392. 


21 .0/22.0 
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Fig. 9.— Changes in the vertical component Bz of the lower boundary field within the region used for the analysis in SectionJ^ 
Panel (a) shows the magnetogram values of Bz in the bin 2 boundary data, for reference. In panels (b)-(e), the changes 
normalized by the provided uncertainties in Bz are displayed for two bin levels for both the magnetofrictional method in panels (b) 
and (c) and for the optimization method in panels (d) and (e). In panels (b) and (d), the contours correspond to the magnetogram 
values of Bz shown in panel (a). 


method, changes are introduced by its version of the prepro¬ 
cessing, which does not impose fixed point-wise constraints. 


and also during the calculation, as described in Section 2.2.1 


Although the Grad-Rubin methods preserve the values of 
the vertical field in the boundary, they do introduce changes 
in the horizontal field in the initial calculation of boundary 
values of the force-free parameter a. Changes are also in¬ 
troduced during solution of the boundary value problem (see 
Section 2.23\ . The boundary values of the horizontal field are 
ignored a priori over one polarity of the vertical field. 


Figures [9| and [T0| show the spatial distribution of the nor¬ 
malized changes and Ah across the lower boundary. In 
Figure the changes in the vertical component Bz of the 
vector magnetogram boundary data introduced by the mag¬ 
netofrictional and optimization methods are shown for both 
the bin 2 and bin 8 calculations. The figure indicates that 
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the magnetofrictional method changes Bz by similar amounts 
(on a normalized basis) for both resolutions, with the more 
highly resolved boundary changed slightly less. The sign of 
Az correlates with the polarity, and indicates that the mag¬ 
netofrictional code systematically reduces the magnitude of 
Bz across the full field of view. The optimization method, in 
contrast, is seen to change the more highly resolved (bin 2 ) 
case significantly less than the lower-resolution (bin 8 ) case, 
on a normalized basis, without any apparent dependence on 
polarity. 


Figure 10 illustrates the changes introduced in the mag¬ 
nitude of the horizontal component Bh of the vector mag¬ 
netogram data, for all methods when using the bin 2 bound¬ 
ary data. This figure shows that A^ is about the same for all 
methods, with the exception for the magnetofrictional method 
where it is lower by about 25% due to the limits imposed dur¬ 
ing the preprocessing step. The normalized changes in Bh 


























Fig. 10.— Changes in the magnitude Bh of the horizontal components of the lower boundary field, for all bin 2 calculations. 
Panels (a) and (b) show the normalized changes for the optimization and magnetofrictional methods. The remaining panels show 
both the P and N solutions for all three Grad-Rubin codes, with panels (c) and (d) for CFIT, panels (e) and (f) for XTRAPOL, and 
panels (g) and (h) for FEMQ. All changes are normalized by the provided uncertainties in the horizontal components of B. Contours 
of Bz are overplotted for reference. 


19 



























are more prominent in the lower portion of the field of view 
where both the values of Bh are the weakest and the uncer¬ 
tainties are the highest. For the Grad-Rubin methods, while 
there are some differences between the P and N solutions in 
the locations of the more prominent changes to B/j, there is a 
good general correspondence between the different methods, 
for a given polarity. This is expected on the basis of the Grad- 
Rubin method: if the different solutions (for a given polarity) 
have a similar connectivity, then field lines carrying strong 
currents will reconnect to the boundary in the opposite polar¬ 
ity at similar locations, and introduce the largest changes in 
the horizontal field at these specific locations. 

4. Discussion and Conclusions 


and computer hardware due to the size (in pixels) of the mod¬ 
eling domain and the corresponding time to completion, and 
thus we report only on the results for the boundary data at bin 
levels 2 to 16. In total, 71 different solution data cubes are 
examined. 

Considered as a group, the physical quantities of interest 
of the NLFFF solutions span a wide range. As calculated 
within the common volume V chosen for analysis, estimated 
total energies E range from 1.08-1.50x 10^® J, estimated free 
energies Ef range from 2% to 24% of the potential field en¬ 
ergies, and estimated relative magnetic helicities Hm range 
from -0.92x 10^® Wb^ to 5.2x 10^® Wb^. The broad ranges in 
these estimates of physical quantities are affected by several 
factors. The three most significant factors are as follows: 


Nonlinear force-free fields (NLFFFs) are increasingly used 
to model the magnetic structure in the Sun’s corona. An ear¬ 
lier series of studies (e.g., Schrijver et al. 12006) Metcalf et al. 


2008 Schrijver et al.|2008 De Rosa et al.||20^ Wheatland 

& Gilchristji2013| l showed that the NLFFF solution methods 

regularly achieved success on fields with known solutions and 
with appropriate boundary conditions, but still encountered 
difficulties in their application of the model to solar data. In 


this paper we re-examine one aspect identified in De Rosa 
et al. ( |2009| l believed to affect the reliability of NLFFF mod¬ 
eling: the influence of the spatial resolution of the boundary 
vector magnetogram data on the NLFFF solutions. 


The boundary data used for this experiment are a sequence 
of vector magnetograms with different spatial resolutions con¬ 
structed from a single normal-map scan in December 2007 of 
AR 10978 by the Hinode/SOT-SR The vector magnetograms 
are produced by rebinning the observed polarization spec¬ 
tra by factors ranging from 2 to 16, performing a Milne- 
Eddington inversion on these rebinned spectra, and then re¬ 
solving the 180° ambiguity in each case. This processing re¬ 
sults in vector magnetograms having a spatial resolution rang¬ 
ing from O'.'SO to 4767, a range of approximately one half to 
one sixteenth of the native Hinode/SOT-SP normal-map res¬ 
olution. The procedure of rebinning the polarization spec¬ 
tra (rather than the inverted vector magnetogram data) is in¬ 
tended to mimic observations by spectrographs with different 
intrinsic resolutions. AR 10978 was a relatively small and 
isolated region at the time of the observation, and as a result 
was selected for this experiment because of the presumption 
that much of the current-carrying field lines lie in the volume 
above the Hinode field of view. The region exhibited some 
magnetic complexity and was flare productive. 

Five different codes implementing solution of the NLFFF 
model with three different methods (the optimization, the 
magnetofrictional, and Grad-Rubin methods) are applied to 
the set of vector magnetogram data produced from the re¬ 
binned polarization spectra. Solutions are calculated using 
these data, spanning the set of spatial resolutions and (for the 
Grad-Rubin methods) for the two choices of polarity in the 
boundary conditions (P and N) on the force-free parameter a. 
Although the Hinode/SOT-SP polarization spectra were also 
processed at the intrinsic instrumental resolution (i.e., with¬ 
out any rebinning), these data proved problematic for codes 


1. The input vector magnetogram data exhibit differences 
across the resolution levels, as shown in Section |2T| 
As a result, some variation in the subsequent NLFFF 
extrapolations, the lower bounds of which are con¬ 
strained by the vector data, is expected. The most gen¬ 
eral trend in the results is that as the spatial resolu¬ 
tion of the boundary data increases, the NLFFF met¬ 
rics improve. In particular, using more highly resolved 
boundary data usually results in NLFFF solutions that 
are both more force- and divergence-free, as found by 
calculating the metrics (CWsin0), (|/i|), and as 
well as by performing the Helmholtz decomposition 


discussed in Section 3.1 Additionally, improvement in 
the NLFFF metrics is correlated with larger Ef in most 
instances. Measurements of Hm, on the other hand, 
do not show any discernible trends with spatial resolu¬ 
tion. Although values of Hm from most methods agree 
to within a factor of two (excepting the optimization 
method, for which values of Hm are clustered around 
zero), the lack of any trend with resolution suggests that 
Hm is difficult to determine from NLFFF solutions. 

2. Vector magnetogram data determined from photo- 
spheric spectral lines are not force-free, and thus are 
not immediately applicable to NLFFF modeling. As 
a result, the NLFFF extrapolation codes necessarily 
change the provided vector boundary data to be more 
compatible with the NLFFF model, resulting in bound¬ 
ary data that are in a more force-free state, much as 
what is expected at the base of the corona. The mea¬ 
surement uncertainties provided with the data are used 
to guide how these changes at the boundary occur, and 
changes to the field (particularly the horizontal compo¬ 
nents) are substantial when compared with the speci¬ 
fied uncertainties, even when higher spatial resolution 
boundary data are used (though the specified uncertain¬ 
ties may be underestimates of the true uncertainties), as 


discussed in Section 3.2 Different extrapolation codes 
implement the changes in different ways, resulting in 
solutions from the different codes, while qualitatively 
similar in appearance (e.g., as in the field-line render¬ 
ings of Fig. 1^, that have quantitatively different char¬ 
acteristics (e.g., the energy estimates shown in Fig. [^. 
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3. Variations from method to method in the cases pre¬ 
sented here may be as large as the effects of spatial 
resolution. For example, the optimization and mag- 
netofrictional methods typically contain about 2 to 5 
times more free energy than the group of Grad-Rubin 
methods. These variations result from a combination 
of effects, including not only the different treatments of 
the vector magnetogram data (point above), but also 
the different conditions imposed on the lateral and top 
boundaries, how unbalanced flux and currents are han¬ 
dled, the size of the held of view relative to the Hux and 
current systems important for NLFFF modeling, and 
the presence or absence of solenoidal errors. Variations 
amongst methods seen in earlier NLFFF extrapolation 
studies (e.g., Schrijver et al.|2006 Metcalf et al.|200^ 


[Schrijver et al.||2008[ De Rosa et al.||2009l were 

larly wide-ranging. 


We conclude that most NLFFF solutions appear more con¬ 
sistent with both the force- and divergence-free conditions 
when more highly resolved vector magnetogram data are 
used. Higher resolution boundary data are associated with 
NLFFF solutions having larger amounts of free energy. How¬ 
ever, the use of more highly resolved boundary data does not 
by itself guarantee a more internally consistent solution. From 
a pragmatic perspective, given the spreads in E, Ef, and Hm 
from the series of NLFFF models shown here and past experi¬ 
ence with constructing NLFFF solutions, we recommend that 
users perform the following checks before a NLFFF model is 
employed in a scientihc setting: 

1. Check metrics such as (CWsin 6), {\fi\), and ^ and as¬ 
sess the degree to which a held is force- and divergence- 
free. In addition, performing a Helmholtz decomposi¬ 
tion on a NLFFF model seems to be a useful way to de¬ 
termine the degree to which the total and free energies 
in the model arise from residual errors in the divergence 
of B. These errors may be signihcant and may call into 
question the accuracies of the free energy estimates, as 
shown in Section 13.1.21 

2. Additionally, either verify that there is good agreement 
between the modeled held lines and the resulting EUV 
and X-ray loop trajectories before using any estimates 
of physical quantities from NLFFF solutions can be re¬ 
lied upon, or use the observed coronal loops to place 
additional constraints during the NLFFF solution pro¬ 
cess. 
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A. Energy Decomposition Details 


This appendix provides more details related to the Helmholtz decomposition of the magnetic energies discussed in Section 3.1.2 
Values of all comTOnents in the decomposition of the magnetic energy E, as defined in Equation (|^, for all of the NLFFF solutions 


a]TOn 

are listed in TableW The tilde over an energy indicates normalization with respect to the total energy, so that, e.g., i?o,s = Eq^s/E. 

There is a clear distinction in results obtained with the magnetofrictional and optimization methods, which seek to minimize 
departures from solenoidality during computation, and the Grad-Rubin implementations, which explicitly solve for a solenoidal 
(divergence-free) field. In the following discussion, we compare the non-solenoidal contributions to the total energy to the Ej ^ 
component, because Ej s is equivalent to the free energy Ej in a perfectly solenoidal field and it is Ef that holds significant physical 
interest. 


The Grad-Rubin code solutions have non-solenoidal contributions which are, in most cases, at least one order of magnitude 
smaller than The Grad-Rubin codes also show little difference in results between the P and N solutions. For the CFIT 
solutions, the calculations using the bin 4 boundary data exhibit the largest non-solenoidal contributions among the CFIT results for 
different resolutions, with values of i?j_ns and | E^i^ \ that lie above the trend established by the CFIT solutions for the other resolution 
levels. Excepting this case, the non-solenoidal errors decrease with increasing resolution from 21% (for bin 16 boundary data) to 1% 
(for bin 2 data) of the free energy Ej^^. The solutions obtained with the XTRAPOL code have dominant non-solenoidal contributions 
from |i?mix| that decrease with increasing resolution from 33% (for bin 16 data) to 1% (for bin 2 data) of Ej^^- The solutions obtained 
with the FEMQ code lie, on average, between the CFIT and the XTRAPOL results, except for the cases using the bin 4 and bin 6 
boundary data. The FEMQ solution for bin 6 boundary data has non-solenoidal contributions of i?j,ns/£'j,s = 3 x 10“^, which are 
the smallest amongst all methods and bin levels. These contributions are also smaller than the non-solenoidal contribution from their 
corresponding potential fields. 


The magnetofrictional solutions exhibit non-solenoidal contributions which are a significant fraction of the free energy in about 
half of the cases. For example, the case using the bin 10 boundary data shows a non-solenoidal contribution of about one third of 
i?j s. The non-solenoidal contributions decrease with increasing resolution down to 6% of the free energy for the case where bin 2 
data are used. The magnetofrictional solutions are affected by non-solenoidal contributions with energies a significant fraction of 
the nominal free energy, and, on average, one order of magnitude larger than for the Grad-Rubin codes. The optimization method 
solutions exhibit the largest non-solenoidal contributions, and the |£itnix| term is larger than the nominal free energy i?j s at most bin 
levels. The mixed term generally increases in size with resolution, from a factor of 1.16 (for bin 16 data) to a factor of 1.52 (for bin 2 
data) of the nominal free energy Ej^s- 


The solenoidal properties of the evolutionary methods are found to improve when the relaxation parameters are adjusted. For 
instance, by changing the parameter Wd in Equation (4) of Wiegelmann et al. (20121 from 1.0 to 1.5, thereby weighting the divergence 
term more strongly during the minimization process, the solutions obtained by the optimization method are reduced by about an order 
of magnitude, as illustrated in Figure [TT] These solutions, however, are less force-free and have less free energy than the solutions 
analyzed in Section For these divergence-optimized solutions, the metrics (CW sin 9) and E/Eq range from 0.46 to 0.52 and 
from 1.02 to 1.08, respectively (cf. Table [^. 

In summary, the Grad-Rubin codes produce solutions with relatively small non-solenoidal contributions, which can be of order 
1% of i?j_s at the highest resolution. The magnetofrictional and optimization methods exhibit significantly larger non-solenoidal 
contributions in energy. Because Equation (8) in |Valori et al. ( 2013| l suggests that non-solenoidal contributions are correlated with 
the magnitude of the current-carrying part of the field, and because the magnetofrictional and optimization method solutions have the 
largest free energies, it is perhaps not surprising that they are characterized by the largest non-solenoidal errors. Irrespective of this 
effect, further reduction of the non-solenoidal errors seems likely to increase the reliability of free energy estimates from solutions 
obtained with the optimization and magnetofrictional methods. 


Finally, we note that the non-solenoidal mixed term i^mix is negative in a number of the NLFFF solutions and may partially cancel 
Fij ns in Equation Large, negative contributions from E^^ can in principle lead to non-physical solutions, with negative free 
energy. This may occur for methods that do not explicitly impose V • B = 0, and especially if non-preprocessed magnetograms 
are used. This was the case for some of the solutions obtained with the magnetofrictional and optimization methods presented in 
Schrijver et al. (2008 1 (cf. Table 1 in that paper). In Valori et al. (2013[) it is shown, for one particular case, how the large negative 


values of E„,ix and non-physical solutions arise (mostly) because of the inconsistency of the vector magnetogram data with the 
force-free equations. 
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Fig. 11.— The solenoidal decomposition for the optimization solution where reducing divergence errors is weighted more strongly 
than reducing the Lorentz force. 


23 




Table 4 

NLFFF Magnetic Energy Decomposition 


Method/Code 

Bin 

i^o,s(xlO-') 

Sj.s(xl 0 - 2 ) 

£^o,ns(xlO^) 

£;j,„s(xl0-3) 

i^mix (xlO 2) 

Optimization 

16 

8.90 

4.93 

4.73 

2.94 

5.75 


14 

8.78 

5.88 

4.15 

3.94 

5.88 


12 

8.48 

5.38 

3.37 

6.03 

9.21 


10 

8.71 

5.78 

3.01 

4.69 

6.63 


8 

8.50 

6.06 

1.98 

4.76 

8.49 


6 

8.47 

7.61 

1.30 

5.16 

7.13 


4 

8.29 

6.69 

0.76 

6.66 

9.75 


3 

8.34 

8.16 

0.50 

6.94 

7.70 


2 

8.08 

7.25 

0.30 

9.07 

11.0 

Magnetofrictional 

16 

8.65 

11.5 

8.17 

16.3 

0.22 


14 

8.68 

12.0 

6.42 

17.9 

-0.65 


12 

8.55 

11.6 

6.53 

21.9 

0.68 


10 

8.37 

12.1 

4.64 

17.8 

2.33 


8 

8.79 

10.5 

3.25 

36.5 

- 2.12 


6 

9.16 

9.3 

2.53 

20.9 

-3.02 


4 

8.89 

11.0 

1.32 

26.8 

-2.58 


2 

9.11 

8.6 

0.43 

5.1 

-0.26 

CFIT (P / N) 

16 

9.60 / 9.67 

4.77/4.15 

15.5 /15.5 

1.19/1.18 

-1.02/-1.07 


14 

9.47 / 9.66 

5.97/4.12 

11.0 /ll.l 

0.91/0.90 

-0.88/-0.90 


12 

9.54/9.57 

5.07/4.72 

9.77/ 9.84 

0.76 / 0.74 

-0.61 /-0.64 


10 

9.49/9.51 

5.52/5.26 

7.56/ 7.59 

0.56/0.58 

-0.53/-0.51 


8 

9.48/9.51 

5.48/5.16 

5.28/ 5.26 

0.38/0.39 

-0.39/-0.36 


6 

9.37/9.51 

6.44 / 5.07 

3.62/ 3.68 

0.27 / 0.27 

-0.25/-0.24 


4 

9.04 / 9.09 

7.56/6.33 

1.84/ 1.88 

1.45/1.93 

1.84/ 2.55 


3 

9.46 / 9.45 

5.50/5.53 

1 . 20 / 1.20 

0.09 / 0.09 

-0.09/-0.09 


2 

9.40/9.51 

6.02/4.95 

0.56/ 0.58 

0.04 / 0.05 

-0.05/-0.05 

XTRAPOL (P / N) 

16 

9.78/9.76 

3.07/3.31 

9.47/9.14 

4.15/3.80 

-1.01 /-1.05 


14 

9.75 / 9.72 

3.29 / 3.60 

6.92/6.71 

3.34/3.13 

-0.890/-0.936 


12 

9.68/9.73 

3.67/3.21 

6.01/5.83 

2.51/2.28 

-0.603/-0.615 


10 

9.59/9.66 

4.60/3.88 

4.76/4.69 

2.02/1.89 

-0.560/-0.575 


8 

9.58/9.63 

4.59/4.07 

3.45 / 3.39 

1.37/1.28 

-0.393 /-0.389 


6 

9.55/9.64 

4.72 / 3.84 

2.42 / 2.40 

0.95 / 0.90 

-0.281 /-0.275 


4 

9.48 / 9.54 

5.28/4.77 

1.34/1.32 

0.50 / 0.48 

-0.146/-0.149 


3 

9.50/9.59 

5.06/4.22 

0.83/0.82 

0.33/0.31 

- 0.101 /-0.106 


2 

9.50/9.56 

5.08/4.48 

0.40 / 0.39 

0.19/0.19 

-0.048/-0.050 

FEMQ (P / N) 

16 

9.66 / 9.67 

3.35/3.25 

9.01/8.68 

7.97 / 7.49 

-0.117/-0.120 


14 

9.73/9.68 

3.08/3.53 

6.92/6.71 

4.39/4.15 

-0.454/-0.461 


12 

9.67 / 9.67 

3.45 / 3.41 

5.91/5.71 

3.92/3.67 

-0.205/-0.197 


10 

9.59/9.63 

4.24 / 3.79 

4.75/4.66 

3.20/3.05 

-0.183/-0.199 


8 

9.55/9.63 

4.68 / 3.94 

3.45 / 3.39 

1.76/1.63 

-0.272 / -0.275 


6 

9.53/9.62 

4.64 / 3.75 

2.44/2.41 

2.07 / 2.04 

0.001 /- 0.001 


4 

9.48 / 9.57 

5.18/4.32 

1.33/1.32 

9.48 / 9.26 

0.003 / 0.001 


3 

9.49 / 9.58 

5.11/4.28 

0.83/0.82 

3.91/3.72 

-0.064/-0.066 


2 

9.45 / 9.54 

5.51/4.57 

0.40 / 0.39 

4.45/4.38 

0.023 / 0.024 


24 



REFERENCES 

Aly, J. J. 1984, ApJ, 283, 349 

—. 1989, Sol. Phys., 120, 19 

Amari, T., & Aly, J.-J. 2010, A&A, 522, A52 

Amari, T., Boulmezaoud, T. Z., & Aly, J. J. 2006, A&A, 446, 
691 

Boulmezaoud, T. Z., & Amari, T. 2000, Z. angew. Math. 
Phys., 51,942 

Chandra, R., Dagum, L., Kohr, D., Maydan, D., McDonald, 
J., & Menon, R. 2001, Parallel Programming in OpenMP 
(San Francisco, CA, USA: Morgan Kaufmann Publishers 
Inc.) 

Chodura, R., & Schlueter, A. 1981, Journal of Computational 
Physics, 41, 68 

Crouch, A. D. 2013, Sol. Phys., 282, 107 

De Rosa, M. L., Schrijver, C. J., Barnes, G., Leka, K. D., 
kites, B. W., Aschwanden, M. J., Amari, T., Canou, A., 
McTiernan, J. M., Regnier, S., Thalmann, J. K., Valori, 
G., Wheatland, M. S., Wiegelmann, T., Cheung, M. C. M., 
Conlon, P. A., Fuhrmann, M., Inhester, B., & Tadesse, T. 
2009, ApJ, 696, 1780 

del Toro Iniesta, J. C. 2003, Introduction to Spectropolarime- 
try (Cambridge, UK: Cambridge University Press) 

Fuhrmann, M., Seehafer, N., & Valori, G. 2007, A&A, 476, 
349 

Fuhrmann, M., Seehafer, N., Valori, G., & Wiegelmann, T. 
2011, A&A, 526, A70 

Georgoulis, M. K. 2012, Sol. Phys., 276, 423 

Grad, H., & Rubin, H. 1958, in Peaceful Uses of Atomic En¬ 
ergy: Proc. Second UN International Atomic Energy Con¬ 
ference, Vol. 31 (Geneva: UN), 190 

Gropp, W., Lusk, E., & Skjellum, A. 1999, Using MPI: 
Portable Parallel Programming with the Message Passing 
Interface, 2nd edn. (Cambridge, MA: MIT Press) 

Henney, C. J., Keller, C. U., & Harvey, J. W. 2006, in Astro¬ 
nomical Society of the Pacific Conference Series, Vol. 358, 
Astronomical Society of the Pacific Conference Series, ed. 
R. Casini & B. W. kites, 92 

Henney, C. J., Keller, C. U., Harvey, J. W., Georgoulis, M. K., 
Hadder, N. L., Norton, A. A., Raouafi, N.-E., & Toussaint, 
R. M. 2009, in Astronomical Society of the Pacific Con¬ 
ference Series, Vol. 405, Solar Polarization 5: In Honor 
of Jan Stenflo, ed. S. V. Berdyugina, K. N. Nagendra, & 
R. Ramelli, 47 


Hoeksema, J. T., Liu, Y., Hayashi, K., Sun, X., Schou, J., 
Couvidat, S., Norton, A., Bobra, M., Centeno, R., Leka, 
K. D., Barnes, G., & Turmon, M. 2014, Sol. Phys., 289, 
3483 

Jackson, J. D. 1999, Classical Electrodynamics, 3rd edn. 
(New York, NY: Wiley) 

Leka, K. D., & Barnes, G. 2012, Sol. Phys., 277, 89 

Leka, K. D., Barnes, G., Crouch, A. D., Metcalf, T. R., Gary, 
G. A., Jing, J., & Liu, Y. 2009, Sol. Phys., 260, 83 

Leka, K. D., Barnes, G., Gary, G. A., Crouch, A. D., & Liu, 
Y. 2012, Sol. Phys., 276, 441 

kites, B. W., Akin, D. L., Card, G., Cruz, T., Duncan, D. W., 
Edwards, C. G., Elmore, D. E, Hoffmann, C., Katsukawa, 
Y, Katz, N., Kubo, M., Ichimoto, K., Shimizu, T., Shine, 

R. A., Streander, K. V, Suematsu, A., Tarbell, T. D., Title, 
A. M., & Tsuneta, S. 2013, Sol. Phys., 283, 579 

Malanushenko, A., Schrijver, C. J., DeRosa, M. L., & Wheat- 
land, M. S. 2014, ApJ, 783, 102 

McClymont, A. N., Jiao, L., & Mikic. 1997, Sol. Phys., 174, 
191 

Metcalf, T. R., De Rosa, M. L., Schrijver, C. J., Barnes, 
G., van Ballegooijen, A. A., Wiegelmann, T., Wheatland, 
M. S., Valori, G., & McTtiernan, J. M. 2008, Sol. Phys., 
247, 269 

Metcalf, T. R., Jiao, L., McClymont, A. N., Canfield, R. C., 
& Uitenbroek, H. 1995, ApJ, 439, 474 

Metcalf, T. R., Leka, K. D., Barnes, G., kites, B. W., Geor¬ 
goulis, M. K., Pevtsov, A. A., Balasubramaniam, K. S., 
Gary, G. A., Jing, J., Li, J., Liu, Y, Wang, H. N., 
Abramenko, V., Yurchyshyn, V, & Moon, Y.-J. 2006, 
Sol. Phys., 237, 267 

Mikic, Z., & McClymont, A. N. 1994, in Astronomical Soci¬ 
ety of the Pacific Conference Series, Vol. 68, Solar Active 
Region Evolution: Comparing Models with Observations, 
ed. K. S. Balasubramaniam & G. W. Simon, 225 

Molodenskii, M. M. 1969, SvA, 12, 585 

Parker, E. N. 1996, ApJ, 471, 485 

Regnier, S. 2013, Sol. Phys., 288, 481 

Roumeliotis, G. 1996, ApJ, 473, 1095 

Scherrer, P. H., Bogart, R. S., Bush, R. L, Hoeksema, J. T., 
Kosovichev, A. G., Schou, J., Rosenberg, W., Springer, L., 
Tarbell, T. D., Title, A., Wolfson, C. J., Zayer, L, & MDI 
Engineering Team. 1995, Sol. Phys., 162, 129 

Schou, J., Scherrer, P. H., Bush, R. L, Wachter, R., Couvidat, 

S. , Rabello-Soares, M. C., Bogart, R. S., Hoeksema, J. T., 
Liu, Y, Duvall, T. L., Akin, D. J., Allard, B. A., Miles, 
J. W., Rairden, R., Shine, R. A., Tarbell, T. D., Title, A. M., 
Wolfson, C. J., Elmore, D. E, Norton, A. A., & Tomczyk, 
S. 2012, Sol. Phys., 275, 229 


25 



Schrijver, C. J., De Rosa, M. L., Metcalf, T., Barnes, G., Lites, 
B., Tarbell, T., McTiernan, J., Valori, G., Wiegelmann, T., 
Wheatland, M. S., Amari, T., Aulanier, G., Demoulin, R, 
Fuhrmann, M., Kusano, K., Regnier, S., & Thalmann, J. K. 
2008, ApJ, 675, 1637 

Schrijver, C. J., De Rosa, M. L., Metcalf, T. R., Liu, Y., Mc¬ 
Tiernan, J., Regnier, S., Valori, G., Wheatland, M. S., & 
Wiegelmann, T. 2006, Sol. Rhys., 235, 161 

Semel, M., & Skumanich, A. 1998, A&A, 331, 383 

Skumanich, A., & Lites, B. W. 1987, ApJ, 322, 473 

Thalmann, J. K., Tiwari, S. K., & Wiegelmann, T. 2013, ApJ, 
769, 59 

Thalmann, J. K., & Wiegelmann, T. 2008, A&A, 484, 495 

Tsuneta, S., Ichimoto, K., Katsukawa, Y, Nagata, S., Otsubo, 
M., Shimizu, T., Suematsu, Y, Nakagiri, M., Noguchi, 
M., Tarbell, T., Title, A., Shine, R., Rosenberg, W., Hoff¬ 
mann, C., Jurcevich, B., Kushner, G., Levay, M., Lites, 
B., Elmore, D., Matsushita, T., Kawaguchi, N., Saito, H., 
Mikami, L, Hill, L. D., & Owens, J. K. 2008, Sol. Rhys., 
249, 167 

Valori, G., Demoulin, R, & Rariat, E. 2012a, Sol. Rhys., 278, 
347 

Valori, G., Demoulin, R, Rariat, E., & Masson, S. 2013, A&A, 
553, A38 

Valori, G., Green, L. M., Demoulin, R, Vargas Dominguez, 
S., van Driel-Gesztelyi, L., Wallace, A., Baker, D., & 
Euhrmann, M. 2012b, Sol. Rhys., 278, 73 

Valori, G., Kliem, B., & Euhrmann, M. 2007, Sol. Rhys., 245, 
263 

Valori, G., Kliem, B., Torok, T., & Titov, V. S. 2010, A&A, 
519, A44 

Wheatland, M. S. 2007, Sol. Rhys., 245, 251 

Wheatland, M. S., & Gilchrist, S. A. 2013, J. Rhys. Conf. Ser., 
440, 012037 

Wheatland, M. S., Sturrock, R. A., & Roumeliotis, G. 2000, 
ApJ, 540, 1150 

Wiegelmann, T., & Inhester, B. 2010, A&A, 516, A107 

Wiegelmann, T., Inhester, B., & Sakurai, T. 2006, Sol. Rhys., 
233,215 

Wiegelmann, T, & Sakurai, T. 2012, Liv. Rev. Solar Rhys., 9, 
5 

Wiegelmann, T., Thalmann, J. K., Inhester, B., Tadesse, T., 
Sun, X., & Hoeksema, J. T. 2012, Sol. Rhys., 281, 37 


This 2-column preprint was prepared with the AAS lATgX macros v5.2. 



